
% !Mode:: "TeX:UTF-8" 

\BiChapter{面向PDOP定位构型的构型控制方法}{Configuration control for PDOP coverage in cooperative localization}

\BiSection{引言}{Introduction}

集群中个体对自身状态的感知是集群在许多实际任务中实现自主功能的前提和基础。
%其中位置状态在许多集群任务中起到了至关重要的作用。
当前集群系统对位置状态的感知和获取途径主要是依赖于GNSS信息。然而在复杂环境中，GNSS定位服务系统的可靠性会严重降低，甚至不可使用。一方面，卫星导航信号经过远距离的空间衰减，到达地面用户时强度低，穿透性差，在隧道、室内、深海等典型环境中无法使用。另一方面，在对抗环境中，GNSS信号容易受到干扰、屏蔽、诱导和欺骗。因此为GNSS拒止环境提供可替代的定位服务具有重要的研究价值。

协同定位由于主要依赖多智能体内部由协作产生的相对观测数据，具有较高的可靠性，因此协同定位方法为多智能体系统在拒止环境中的定位任务提供了可靠的替代解决方案。虽然已经有大量的研究结果表明，协同定位的估计误差精度、以及误差的增长速度与多智能体系统内部的相对观测构型有关。但是，针对应当如何优化系统构型从而提升定位性能的问题尚且缺乏较为系统的研究。文献\cite{1986Optimal,2000GPS}针对GNSS系统下卫星与地面用户的相对构型与GNSS定位精度之间的关系进行了分析。文献\cite{20113D,2008Optimal,1999Dilution}则探讨了伪卫星的位置对增强GNSS定位几何构型的方法。然而相关研究成果存在如下几点局限性：

(1) 锚节点保持预定义的周期性运动(卫星)，或者保持静止(地基伪卫星基站)，因此定位性能无法实时在线调整。

(2) 得出的优化构型往往是针对特定数量锚节点的特解。无法扩展至任意数量的锚节点集群。

(3) 构型设计和分析的过程往往是集中式的。

基于上述问题，本章的主要目的是定量地研究构型控制与协同定位之间的耦合关系，即通过设计构型控制方法，使得移动协同定位系统具备实时调整定位性能、拓展定位服务范围、自适应任意系统规模等特点。本章使用GNSS研究领域中可以反映定位误差受构型影响的位置精度影响因子(PDOP)作为构型控制指标，设计了一种旨在优化GNSS拒止区域内PDOP场覆盖强度的分布式构型控制方法。基于覆盖控制框架\citeup{cortes2004coverage,arslan2019statistical,zhong2011distributed,cassandras2005sensor}，将面向PDOP定位构型的协同控制问题转化为了PDOP场在拒止区域的覆盖控制问题，通过严格的理论推导得出了覆盖指标对每一个多智能体节点位置状态的导数表达式，并以此为基础设计了基于梯度下降更新的协同定位构型分布式构型控制方法。最后通过数值仿真验证了该构型控制方法对大范围GNSS拒止区域中静态热点定位、节点丢失和加入、切换热点等三类定位场景下构型变换对协同定位过程的提升效果，验证了所提方法在理论上的正确性，以及在实际应用中对系统初始构型的鲁棒性。
%一种常见方法是使用伪基站对GNSS信号进行增强、中继和补充。然而在当前对定位服务系统的研究中，地基伪基站总是保持静止，空基伪基站则总是当前基于伪基站的移动定位服务系统（空基和地基）


%对GNSS受限拒止区域的定位导航问题，已有的研究工作目前主要聚焦于以下两类方法：
%\begin{itemize}
%	\item[1):] 在受限拒止区域部署伪基站。伪基站是一类可以接收GNSS信号，并向其周边区域用户发射增强后的类GNSS信号的设备。由于伪基站一般部署于地面或者近地空间，不受电离层干扰影响，受信号在自由空间衰减的影响较少，因此伪基站周边区域导航信号的强度比卫星的导航信号强度将高出多个数量级，从而有效地克服信号视线遮挡，为存在较多障碍物的环境提供更稳定的导航服务。在GNSS信号完全拒止的环境里，可以通过人工校准获取基站的精准位置。此类方法的典型应用为LOCATA地基定位系统。该系统在美国空军的白沙导弹靶场测试中实现了误差仅为18 $cm$的定位精度，其理论精确可达5 $cm$。
%	\item[2):] 配备新型传感器以获取额外的导航定位信息，从而通过融合不同定位原理的位置信息源，实现稳定的导航定位功能。此类技术中涉及的信息源包括距离信息、信号频率、信号强度、信号指纹、线性/角度加速度信息等。由此使用到的传感器包括射频类（如WLAN、WIFI、蓝牙、UWB以及RFID等），超声波，红外射线，惯性测量单元，视觉辅助系统以及可见光通信等。针对不同类型的信息源，其对应的定位原理各不相同，由此衍生的定位技术包括三角定位法、三边测量法、数据匹配方法等。然而仅依靠少量的测量数据无法获得稳定的导航定位效果，因此基于优化方法的融合定位（如融合滤波，图优化等）逐渐成为首选。此类方法的代表为同时定位与建图（SLAM）以及协同导航（CL）等热点研究。
%\end{itemize}

\BiSection{问题描述}{Problem formulation}
\BiSubsection{场景与数学符号定义}{Scenario and notations}
如图~\ref{fig. PDOP_notations}所示，任务空间被定义为一处由二维凸平面$\vect{\Omega} \in \mathbb{R}^2$通过垂直投影形成的三维有界空间$\vect{F} \in \mathbb{R}^3$。 假设由于信号干扰或者环境限制，存在于该任务空间地表的各类单位均无法获取GNSS导航定位信息，本章研究如何利用一组空基移动节点为任务空间中存在的地面用户提供定位服务。与卫星导航系统类似，本章假设空基节点的自定位信息精确已知。空基节点的自定位过程可以通过使用成熟的定位算法，如协同SLAM算法等，从环境中获取必要的定位信息实现自定位功能。然后，空基节点通过携带的导航设备，向地面单位发射存储自身位置信息和时间戳的导航信号。本章假设地表用户均携带合适的信号接收设备，可以正确地接收空基节点发射的信号并解码信号中传输的导航信息。通过测量导航信号的到达时间（Time of Arrival, TOA），地面用户即可正确解算自身位置。

\begin{figure}[htbp]
	\centering
	\includegraphics[width = 0.9\textwidth]{figures_PDOP/NotationDefinition}
	\bicaption[fig. PDOP_notations]{}{任务空间与数学符号}{Fig.$\!$}{The mission space and notation definitions}\vspace{-0.5em}
\end{figure}

任务空间如图~\ref{fig. PDOP_notations}所示，使用三维向量$\vect{S}_i = [s_{ix},s_{iy},s_{iz}]^T$表示第$i$个空基节点在任务空间中的位置。每个节点通过使用导航设备向地表用户提供导航定位服务。假设每个节点有效的导航服务覆盖模型为地表平面上以该节点垂直投影为中心的圆形覆盖区域 $\vect{\Omega}_i$。考虑飞行平台受到的实际物理约束，例如信号发射设备的功率受限等因素，单个节点的有效覆盖范围是受限的（即覆盖半径$r_i$是有界的）且受节点飞行高度的改变而变化。该覆盖模型可以统一表述为如下的数学模型：
\begin{equation} \label{eq. generalCoverageModel}
r_i = r(s_{iz},L)
\end{equation}
式中 $L$ 表征导航设备因受到实际的物理约束而存在的约束变量。例如图~\ref{fig. PDOP_notations}中展示了一类与节点的飞行高度成反比例的覆盖模型，此处的 $L$ 为导航信号由于受到发射功率的限制而能传播的最远距离。作为公式\eqref{eq. generalCoverageModel}的特例，该覆盖模型可表示为，
\begin{equation} \label{eq. reciprocalCoverageModel}
r_i = \sqrt{L^2 - s^2_{iz}}
\end{equation}   
后文的控制策略设计过程并不需要覆盖模型精确的数学表达式，只要求该数学模型对节点的飞行高度连续、可导。本章延用几何学的通用符号集，使用前置符号 $\partial$ 表征一个集合的边界，即采用$\partial \vect{\Omega}_i$ 表示第 $i$ 个节点导航服务有效覆盖区域的地面圆形边界。


\BiSubsection{精度因子}{Dilution of precision}
前一节介绍了协同导航服务的任务空间定义、节点配置以及必要的数学定义。为了评估协同导航系统的服务质量，本节引入卫星导航研究领域的精度因子（Dilution of Precision, DOP）概念。在卫星导航领域，由于卫星一般运行在中高轨道上，对于一个地面用户而言，来自不同卫星的导航信号的伪距测量精度可认为是固定且相同。因此，用户的定位精度可以使用几何精度因子(Geometry DOP, GDOP)来衡量，并且可以使用如下公式计算：
$$GDOP = \sqrt{PDOP^2+TDOP^2}$$
上式中，PDOP是衡量卫星星座和用户所在位置相对构型优劣的位置精度因子，TDOP是衡量用户接收设备与卫星基准时间在时间上偏移误差值的时间精度因子。

由于本章关注多智能体构型对协同定位的定量影响，因此使用位置精度因子PDOP评估移动协同导航系统的服务质量。对于一个地面用户，使用 $\vect{p} = [p_x,p_y]^T$ 代表其在二维平面的位置。定义$\vect{N}^k_{\vect{p}}$为用户$\vect{p}$在$k$时刻观测到的空基节点的集合$\vect{N}^k_{\vect{p}}:=\{ \vect{S}_i | \vect{p} \in \vect{\Omega}_i \}$，使用$n^k_{\vect{p}}$表示集合$\vect{N}^k_{\vect{p}}$中包含的元素数量。则在用户处的PDOP值可以按如下三个步骤进行计算：

\begin{itemize}
	\item[(1)] 将用户的二维坐标扩展到三维空间 $\bar{\vect{p}} = [\vect{p},0]^T$,然后针对其观测到的空基节点集合$\vect{N}^k_{\vect{p}}$ 中的每一个节点 $\vect{S}_i$计算用户到空基节点的单位方向向量：
	\begin{equation} \label{eq. unitDirectionVector}
		\vect{h}_i = \vect{h}_{\vect{p} \rightarrow \vect{S_i}} = \frac{\vect{\bar{p}} - \vect{S}_i }{|| \vect{\bar{p}} - \vect{S}_i ||_2}
	\end{equation} 
	\item[(2)]对$\vect{N}_{\vect{p}}^k$ 中的每一个节点计算其单位方向向量后，将所有的方向向量整合到一个$n^k_{\vect{p}}\times3$的联合观测矩阵中：
	\begin{equation} \label{eq. collectedHMatrix}
		\vect{H} = \left[ \begin{matrix}
		\vect{h}_1^T \\
		\vdots \\
		\vect{h}_{n^k_{\vect{p}}}^T
		\end{matrix} \right]
	\end{equation}
	\item[(3)] 用户$\vect{p}$处空基节点处用户的PDOP可按照如下公式进行计算：
	\begin{equation} \label{eq. calculationPDOP}
		PDOP(\vect{p}, \vect{S}) = tr\left( \left( \vect{H}^T\vect{H} \right)^{-1} \right)
	\end{equation}
	上式中$tr(\cdot)$为求矩阵的迹。
\end{itemize}

由于PDOP本质上是伪距的测量误差到用户自定位误差的缩放因子，因此当移动导航系统的节点个数等硬件配置已经确定后，降低用户所在位置的PDOP指标可以有效地提升协同导航服务的质量,减少用户的定位误差。

与针对单个用户的定位任务有所不同，移动协同定位系统构型控制的目的是为给定热点区域提供主动的定位导航服务，即优化任务空间内所有可能存在的用户接收到的定位服务质量。因此本章待解决的问题是：设计空基节点的构型控制方法，使得任务空间中的目标区域给定后，空基节点具备自动快速机动到目标区域上空的能力，并通过内部交互自动地形成可以提供最优PDOP分布的相对构型。

\BiSection{PDOP场的协同覆盖优化}{PDOP filed coverage formulation}

本节将沿用覆盖控制框架，将前一节提出的面向PDOP指标的构型控制问题，转化为协同覆盖控制框架下的优化问题，以方便后续的控制方法设计。同时还将解决导航设备的相对测量距离有界受限时，PDOP场在地面空间上分布不连续的问题。

为了描述目标区域需要导航服务的地面用户分布，定义一类密度函数$R(\vect{p}):\vect{\Omega} \rightarrow \mathbb{R}$ 来表征和量化任务空间内某一地面坐标需要获取导航服务的重要性。该密度函数在任务空间中是正定且有界的，即$R(\vect{p}) > 0 $ 并且 $\int_{\vect{\Omega}} R(\vect{p}) d \vect{p} < \infty$。一类典型的密度函数如下所示：
\begin{equation} \label{eq. gaussianDensity}
	R(\vect{p}) = exp(-\alpha||\vect{p} - \vect{p}_d||_2)
\end{equation} 
式中$\vect{p}_d$为密度中心点坐标，$\alpha$为密度函数分布的衰减速率。由式\eqref{eq. gaussianDensity}易知，该密度呈高斯分布特征，因此也被称为高斯密度函数。

由此，为了获取任务空间内最优的PDOP覆盖分布，可通过最大化如下的积分指标函数：
\begin{equation} \label{eq. PDOPultility1}
	\mathcal{H}(\vect{S}) = \int_{\vect{\Omega}} R(\vect{p}) f\left( PDOP(\vect{p},\vect{S}) \right) d \vect{p}
\end{equation}
式中$f(\cdot)$是一类单调递减、可微分的连续函数。由于PDOP的数值大小与定位服务质量成反相关，因此该函数的作用是将此关系转变为正相关，以符合将上述指标函数最大化的目的。

然而，在给定空基节点分布后，PDOP场在地面空间上的分布不是连续的。一方面，由于导航设备存在最大通信测量距离约束，地面用户$\vect{p}$在某一时刻所能观测到的空基节点数量$n^k_{\vect{p}}$（即相对测量数量）不是固定不变，而是随着时间变化和空间相对位置的变化而该改变的。由于联合观测矩阵$\vect{H} \in \mathbb{R}^{n^k_{\vect{p}}\times3}$，$\vect{H}^T\vect{H}$为$3\times3$的方阵。但是如果$det(\vect{H}^T\vect{H})<3$，式~\eqref{eq. calculationPDOP}中的取逆操作无效，也因此PDOP值将无法计算。另一方面当空基节点与用户坐标共线或共面时，也将导致$det(\vect{H}^T\vect{H})<3$。 总体而言，当空基节点集群的相对构型出现如下配置时，将会导致PDOP计算过程出现奇异：

\begin{itemize}
	\item[(1)] 用户观测到的节点数目$n^k_{\vect{p}}$小于3，此时$det(\vect{H}^T\vect{H})=0$;
	\item[(2)] 用户观测到的节点个数$n^k_{\vect{p}}\ge3$个，但是所有节点全部与用户坐标共线，此时$det(\vect{H}^T\vect{H})=1$。
	\item[(3)] 用户观测到的节点个数$n^k_{\vect{p}}\ge3$个，但是所有节点都与用户坐标位于共同平面，此时$det(\vect{H}^T\vect{H})=2$;
\end{itemize}

由于PDOP在上述特殊的构型下无定义，因此这些奇异构型会对覆盖指标的优化产生影响。为了避免这些奇异配置，需要将它们与正常配置进行分离。因此定义如下的分段函数：
\begin{equation} \label{eq. modificationGeneralPDOP}
	\begin{aligned}
	& \bar{f}\left( \cdot \right) = 
	\left\lbrace
	\begin{aligned}
	a, \quad	& 	\text{如果\quad} n_{\vect{p}}^k = 1\\
	b, \quad 	&	\text{如果\quad} n_{\vect{p}}^k = 2\\
	c, \quad 	&	\text{如果\quad} n_{\vect{p}}^k \ge 3 \text{\quad 但是\quad} det(\vect{H}^T\vect{H}) < 3 \\ 
	f\left( \cdot \right), \quad	& \text{其他情况}
	\end{aligned} 
	\right. 
	\end{aligned}
\end{equation}
式中$a,b,c$为三个大于零的常数。它们的取值需要远远小于正常节点配置下的PDOP值。虽然奇异配置下的PDOP数值无定义，但是观测到的节点数目越多，对优化其未来时刻的PDOP数值越有利，因此在实际应用中推荐设置$a<b<c$。由于该分段函数在奇异配置下不存在导数，为了方便后续的推导过程，其导数函数被强制设置如下：
\begin{equation} \label{eq. DerivativeModiGeneralPDOP}
	\begin{aligned}
	& \partial \bar{f}\left( \cdot \right) = \\
	& \left\lbrace
	\begin{aligned}
	0, \quad	& 	\text{如果 \quad } n_{\vect{p}}^k = 1,2\\
	0, \quad 	&	\text{如果 \quad} n_{\vect{p}}^k \ge 3 \text{ 但是 \quad} det(\vect{H}^T\vect{H}) < 3 \\ 
	\partial f\left( \cdot \right), \quad	& \text{其他情况}
	\end{aligned} 
	\right. 
	\end{aligned} \quad.
\end{equation}

公式~\eqref{eq. modificationGeneralPDOP}和~\eqref{eq. DerivativeModiGeneralPDOP}对PDOP分布场的修改，使得PDOP场的定义在任务空间内所有地面坐标上有意义。因此，式~\eqref{eq. PDOPultility1}的优化指标函数被替换为
\begin{equation} \label{eq. objectiveFunction}
	\mathcal{H}(\vect{S}) = \int_{\vect{\Omega}} R(\vect{p}) \bar{f} \left( PDOP\left( \vect{p},\vect{S} \right) \right) d\vect{p}
\end{equation}

基于以上的分析，面向PDOP定位构型的构型控制问题可以转化为求解如下的优化问题：
\begin{equation} \label{eq. modifedObejective}
	\begin{aligned}
		\mathop{\max}_{\vect{S}} & \int_{\vect{\Omega}} R(\vect{p}) \bar{f} \left( PDOP\left( \vect{p},\vect{S} \right) \right) d\vect{p} \\
		\text{s.t.} & \qquad \vect{S}_i \in \vect{F} & i = 1,\cdots,N\\
		& \qquad r_i = r(s_{iz}) &i = 1,\cdots,N\\
	\end{aligned}
\end{equation}

\BiSection{分布式控制策略设计}{Distributed control for PDOP configuration}
前两节通过对问题的描述和分析建立了针主动导航服务的优化数学模型，该模型的优化目的是最大化PDOP场在任务空间对地面目标目标用户的覆盖分布。本节将基于梯度下降方法设计一类分布式控制策略。当热点区域通过密度函数给出后，空基节点可以通过该控制策略自动地机动到该热点上方，并为目标区域提供优化PDOP分布。

如果两个节点的有效导航服务的覆盖区域重合，则称此两个节点之间互为临接节点。定义$\vect{N}_{\vect{S}_i}^k:=\{ \vect{S}_j | \vect{\Omega}_i \bigcap \vect{\Omega}_j \neq \varnothing \}$ 为第$i$个节点在第$k$时刻的所有临接节点的集合。在设计构型控制方法之前，首先给出本章的一个合理假设：

\begin{assumption} \label{ass. PDOPneighboorSet}
	在任意时刻，任意节点均可获取其临接节点集合$\vect{N}_{\vect{S}_i}^k$。
\end{assumption}

上述假设可以通过装配必要的硬件设备来实现。例如，如果使用无人机集群作为空基节点载体，则除了前文介绍的必要的导航和自定位设备外，每个无人机还可以携带一款专用的通信设备，用来与临界节点互相传递自身位置信息。实现该假设的一个必要条件为，该专用通信设备的通信距离是导航设备最大服务距离的2倍大小。此时，当两个节点的有效导航覆盖区域相交或相切时，它们彼此之间可以通过专用通信设备建立通信连接，从而各自组建其对应的邻接节点集合。

\BiSubsection{基于梯度下降方法的策略设计}{PDOP configuration control based on gradient-descent method}

受到无线信号网络覆盖控制方法设计的启发，本节将设计一类梯度下降控制器。根据梯度下降法的更新规则，第$i$个节点的位置更新策略为：
\begin{equation} \label{eq. gradienDescentUpdateRule}
	\vect{S}_i^{k+1} = \vect{S}_i^k + \lambda_k \frac{\partial \mathcal{H}(\vect{S})}{\partial \vect{S}_i^k}
\end{equation}
式中$\lambda_k$为更新步长，$\frac{\partial \mathcal{H}(\vect{S})}{\partial \vect{S}_i^k}$为指标函数对于第$i$个节点位置的偏导数。该控制方法的难点在于指标函数对每个节点位置偏导数的计算。因此，后文重点关注于如何计算该偏导数。

\begin{figure}[htbp]
	\centering
	\includegraphics[width = 0.6\textwidth]{figures_PDOP/GradientOnBoundary1}
	\bicaption[fig. PDOP_GradientOnBoundary]{}{指标函数在覆盖边界上的梯度}{Fig.$\!$}{Gradient of the objective on the coverage boundary}\vspace{-0.5em}
\end{figure}

通过观察指标函数在式~\eqref{eq. objectiveFunction}中的定义，
%对于任务空间地面某个坐标$\vect{p}$,计算该点处的PDOP需要用到所有覆盖该坐标的节点位置。因此指标函数中的积分操作无法如文献\cite{cortes2004coverage}分解到单个节点的服务覆盖域中。然而
指标函数对第$i$个节点位置的偏导数却仅在节第$i$个点的有效覆盖区域内存在。这是因为该区域外的地面坐标的PDOP计算不依赖于第$i$个节点的位置。因此指标函数对节点位置的偏微分可以按照如下方式计算：
\begin{equation}
	\frac{\partial \mathcal{H}(\vect{S})}{\partial \vect{S}_i} =
	\frac{\partial}{\partial \vect{S}_i} \left[ \int_{\vect{\Omega}_i}  R(\vect{p}) \bar{f}\left( PDOP(\vect{p},\vect{S}) \right) d \vect{p}  \right]
\end{equation} 
上述偏微分过程由于独特的数学性质限缩了被导函数的积分域，使得被导函数不需要考虑当前节点与其他节点的耦合关系，在一定程度上简化了该问题。然而，被导函数中的积分域$\vect{\Omega}_i$和被积分函数$R(\vect{p}) \bar{f}\left( PDOP(\vect{p},\vect{S}) \right)$均为节点位置$\vect{S}_i$的函数。对此类积分函数的微分过程需要使用对积分函数求导的莱布尼兹法则以及一种“物理拓展”方法\citeup{1973Differentiation}。据此，该偏微分可分解如下：
\begin{equation} \label{eq. overviewPartialDerivative}
	\begin{aligned}
		\frac{\partial \mathcal{H}(\vect{S})}{\partial \vect{S}_i} &= \int_{\vect{\Omega}_i} R(\vect{p}) \frac{\partial \bar{f} \left( PDOP(\vect{p},\vect{S})\right)}{\partial \vect{S}_i} d\vect{p} \\
		& + \int_{\partial \vect{\Omega}_i} R(\vect{p})\bar{f}\left( PDOP(\vect{p},\vect{S}) \right) \vect{v}\cdot\vect{n} dw
	\end{aligned}
\end{equation}
上式将指标函数的偏微分分解成两个部分，其中第一项代表指标函数在第$i$个节点覆盖区域内部的梯度，第二项为指标函数在覆盖区域边界处的梯度。上式右方第二项中$dw$，$\vect{v}$和$\vect{n}$为新引入的变量，它们的“物理”定义可参见图~\ref{fig. PDOP_GradientOnBoundary}。 其中$dw$为覆盖区域边界处积分点坐标$\vect{p}$在边界上的单位弧长，$\vect{n}$为边界在积分点坐标$\vect{p}$处的单位法向量，$\vect{v}$为积分点$\vect{p}$处坐标因节点位置$\vect{S}_i$的变化而导致的积分点坐标在地面的运动速度向量（或矩阵）。

后文将详细介绍指标函数的微分操作分别在覆盖区间内和覆盖区域边界处的计算方法。

\BiSubsection{覆盖区间的梯度计算}{Graident within coverage area}
本节关注式~\eqref{eq. overviewPartialDerivative} 右侧第二项微分的具体计算与实现。由于密度函数$R(\vect{p})$与节点的位置$\vect{S}_i$无关，因此该微分的计算主要在于求解分段函数$\bar{f}(\cdot)$对节点位置的导数关系。根据微分的链式求导法则，
\begin{equation}
	\frac{\partial \bar{f} \left( PDOP(\vect{p},\vect{S})\right)}{\partial \vect{S}_i} = \frac{\partial \bar{f}}{\partial PDOP} \frac{\partial PDOP(\vect{p},\vect{S})}{\partial \vect{S}_i}
\end{equation}

由公式~\eqref{eq. modificationGeneralPDOP}和~\eqref{eq. DerivativeModiGeneralPDOP}的定义可得，分段函数$\bar{f}(\cdot)$对节点位置的导数仅依赖于连续函数$f(\cdot)$的选取。因此上式中第二项$\frac{\partial PDOP(\vect{p},\vect{S})}{\partial \vect{S}_i}$是计算该整个区间内指标函数梯度的关键。然而PDOP在坐标点$\vect{p}$处的计算涉及多次矩阵运算。因此下列三项关于矩阵的求导公式是进一步计算的关键：
\begin{align}
	\frac{\partial}{\partial x} f(\vect{A}(x)) =& tr\left[ \frac{\partial f}{\partial \vect{A}} \frac{\partial \vect{A}}{\partial x} \right] \\
	\frac{\partial }{\partial \vect{A}} tr(\vect{A}) =& \vect{I} \\
	\frac{\partial }{\partial x} \vect{A}^{-1} = & -\vect{A}^{-1} \left( \frac{\partial \vect{A}}{\partial x} \right)  \vect{A}^{-1}
\end{align} 

因此，结合链式法则和上式中给出的矩阵运算求导规则，$\frac{\partial PDOP(\vect{p},\vect{S})}{\partial \vect{S}_i}$的计算可进一步推导如下：
\begin{equation} \label{eq. derivativePDOPSi}
	\begin{aligned}
		& \frac{\partial PDOP(\vect{p},\vect{S})}{\partial \vect{S}_i} = \frac{\partial tr \left[ \left( \vect{H}^T \vect{H} \right)^{-1} \right]}{\partial \vect{S}_i} \\
		& =  tr \left[ \frac{\partial tr \left[ \left( \vect{H}^T \vect{H} \right)^{-1} \right]}{\partial \left( \vect{H}^T \vect{H} \right)^{-1}} \frac{\partial \left( \vect{H}^T \vect{H} \right)^{-1}}{\partial \vect{S}_i}\right] \\
		& =  tr \left[ \vect{I} \frac{\partial \left( \vect{H}^T \vect{H} \right)^{-1}}{\partial \vect{S}_i}\right] \\
		& = tr \left[ -\left( \vect{H}^T \vect{H} \right)^{-1} \frac{\partial \vect{H}^T\vect{H}}{\partial \vect{S}_i} \left( \vect{H}^T \vect{H} \right)^{-1} \right]
	\end{aligned}
\end{equation}

考虑一个地面坐标$\vect{p}\in \vect{\Omega}_i$,由上式可以看出，为了计算在该点处指标函数对节点$\vect{S}_i$的微分，需要获得$\vect{p}$处的联合观测矩阵$\vect{H}$。由式~\eqref{eq. collectedHMatrix}观测矩阵的定义可知，矩阵$\vect{H}$的计算依赖于$\vect{p}$处的观测到的节点集合$\vect{N}_{\vect{p}}^k$。基于假设~\ref{ass. PDOPneighboorSet},集合$\vect{N}_{\vect{S}_i}^k$可以通过假设额外设备获得。分析它们的定义可得如下关系：
\begin{equation}
	\vect{N}_{\vect{p}}^k \subset \vect{N}_{\vect{S}_i}^k \bigcup \{ \vect{S_i} \}
\end{equation}
因此，如果假设~\ref{ass. PDOPneighboorSet}成立，则节点$\vect{S}_i$可以自行计算其其覆盖域内所有坐标点的联合观测矩阵$\vect{H}$.

再次由$\vect{H}$矩阵在式~\eqref{eq. collectedHMatrix}的定义可得：
\begin{equation}
	\begin{aligned}
		\frac{\partial \vect{H}^T\vect{H}}{\partial \vect{S}_i} & = \frac{\partial \vect{h}_1^T\vect{h}_1}{\partial \vect{S}_i} +\frac{\partial \vect{h}_2^T\vect{h}_2}{\partial \vect{S}_i} + \cdots + \frac{\partial \vect{h}_{n_{\vect{p}}^k}^T\vect{h}_{n_{\vect{p}}^k}}{\partial \vect{S}_i} \\
		& = \frac{\partial \vect{h}_i^T\vect{h}_i}{\partial \vect{S}_i}
	\end{aligned}
\end{equation}

带入单位观测向量$\vect{h}_i$：
\begin{equation} \label{eq. derivativehinArea_x}
	\frac{\partial \vect{h}_i^T \vect{h}_i}{\partial s_{ix}} = \frac{1}{||\vect{\bar{p}} - \vect{S}_i||_2^2} \left[  \begin{matrix}
	2(s_{ix} - p_x) 	&	s_{iy} - p_y	&	s_{iz} 	\\
	s_{iy} - p_y		&		0			&		0	\\
	s_{iz}				&		0			&		0	\\
	\end{matrix} \right]
\end{equation}

\begin{equation}\label{eq. derivativehinArea_y}
	\frac{\partial \vect{h}_i^T \vect{h}_i}{\partial s_{iy}} = \frac{1}{||\vect{\bar{p}} - \vect{S}_i||_2^2} \left[  \begin{matrix}
	0				 	&	s_{ix} - p_x	&		0 	\\
	s_{ix} - p_x		&	2(s_{iy} - p_y)	&	s_{iz}	\\
	0					&		s_{iz}		&		0	\\
	\end{matrix} \right]
\end{equation}

\begin{equation} \label{eq. derivativehinArea_z}
	\frac{\partial \vect{h}_i^T \vect{h}_i}{\partial s_{iz}} = \frac{1}{||\vect{\bar{p}} - \vect{S}_i||_2^2} \left[  \begin{matrix}
	0				 	&	0					&	s_{ix} - p_x 	\\
	0					&	0					&	s_{iy} - p_y	\\
	s_{ix} - p_x		&	s_{iy} - p_y		&	2s_{iz}	\\
	\end{matrix} \right]
\end{equation}
其中$\vect{\bar{p}}$为地面而为用户拓展到三维空间的坐标，即$\vect{\bar{p}} = [\vect{p}^T,0]^T$. 

由上述推导过程可总结如下的梯度计算公式：
\begin{equation} \label{eq. computationGradientArea}
	\begin{aligned}
		&\int_{\vect{\Omega}_i} R(\vect{p}) \frac{\partial \bar{f} \left( PDOP(\vect{p},\vect{S})\right)}{\partial \vect{S}_i} d\vect{p} = \\
		&\left[ \begin{matrix}
		\int_{\vect{\Omega}_i} -R(\vect{p}) \frac{\partial \bar{f}}{\partial PDOP} tr \left[ \left( \vect{H}^T \vect{H} \right)^{-1} \frac{\partial \vect{h}_i^T\vect{h}_i}{\partial \vect{S}_{ix}} \left( \vect{H}^T \vect{H} \right)^{-1} \right] d\vect{p} \\ 
		\int_{\vect{\Omega}_i} -R(\vect{p}) \frac{\partial \bar{f}}{\partial PDOP} tr \left[ \left( \vect{H}^T \vect{H} \right)^{-1} \frac{\partial \vect{h}_i^T\vect{h}_i}{\partial \vect{S}_{iy}} \left( \vect{H}^T \vect{H} \right)^{-1} \right] d\vect{p} \\
		\int_{\vect{\Omega}_i} -R(\vect{p}) \frac{\partial \bar{f}}{\partial PDOP} tr \left[ \left( \vect{H}^T \vect{H} \right)^{-1} \frac{\partial \vect{h}_i^T\vect{h}_i}{\partial \vect{S}_{iz}} \left( \vect{H}^T \vect{H} \right)^{-1} \right] d\vect{p}
		\end{matrix} 
		\right]
	\end{aligned}
\end{equation}
式中$\frac{\partial \vect{h}_i^T\vect{h}_i}{\partial \vect{S}_{i}}$在$x,y,z$三个方向的分量可以参照公式~\eqref{eq. derivativehinArea_x}-\eqref{eq. derivativehinArea_z}。

\BiSubsection{覆盖边缘的梯度计算}{Gradient on the boundary}
前一节通过严格的推导过程得出了指标函数在覆盖区域内的梯度项。本节关注式~\eqref{eq. overviewPartialDerivative}中右侧第二项的梯度的计算。如图~\ref{fig. PDOP_GradientOnBoundary}所示，首先将覆盖区域边界处的地面坐标$\vect{p} \in \partial \vect{\Omega}_i$转换到以覆盖区域圆心$[s_{ix},s_{iy}]^T$为原点的极坐标系下，
\begin{align} 
	p_x &= r_icos(\theta) + s_{ix} \label{eq. 2D-3DTransform_x}\\
	p_y &= r_isin(\theta) + s_{iy} \label{eq. 2D-3DTransform_y},
\end{align}
式中$\theta$为相对角度。

根据定义，向量$\vect{n}$是点$\vect{p}$在边界上的垂直法向量。因此可以通过将$\vect{p}$处的切线方向的向量旋转90度来获得。因此首先计算$\vect{p}$处在边界上的切线方向单位向量:
$$\frac{1}{dw}[dp_x,dp_y]^T = \frac{1}{dw} [-r_i\sin(\theta)d\theta,r_i\cos(\theta)d\theta ]^T$$
由此向量$\vect{n}$可以按照下式选取：
\begin{equation}\label{eq. normalVector}
	\vect{n} = \frac{1}{dw} \left[  \begin{matrix} 
	dp_y\\
	-dp_x
	\end{matrix} \right]
	= \frac{d\theta}{dw}\left[ \begin{matrix}
	r_i\cos(\theta) \\
	r_i\sin(\theta)
	\end{matrix} \right] 
\end{equation}

此外，依据向量$\vect{v}$的含义，其本质是代表节点的运动对覆盖边界的影响。因此该变量的获取可以通过求取地面坐标对对单个节点的偏微分，即：
\begin{equation} \label{eq. velocityVector}
	\vect{v} = \left[ \begin{matrix}
	\frac{\partial p_x}{\partial s_{ix}} & \frac{\partial p_y}{\partial s_{ix}} \\
	\frac{\partial p_x}{\partial s_{iy}} & \frac{\partial p_y}{\partial s_{iy}} \\
	\frac{\partial p_x}{\partial s_{iz}} & \frac{\partial p_y}{\partial s_{iz}} 
	\end{matrix}
	\right] = \left[ \begin{matrix}
	1	& 	0\\
	0	&	1\\
	\frac{\partial r}{\partial s_{iz}}\cos(\theta) 	& \frac{\partial r}{\partial s_{iz}}\sin(\theta)
	\end{matrix}
	\right]
\end{equation}

将式~\eqref{eq. velocityVector}和式~\eqref{eq. normalVector}带入式~\eqref{eq. overviewPartialDerivative}，则指标函数在覆盖边界处的梯度可按照如下公式进行计算：
\begin{equation} \label{eq. computationGradientBoundary}
	\begin{aligned}
		&\int_{\partial \vect{\Omega}_i} R(\vect{p})\bar{f}\left( PDOP(\vect{p},\vect{S}) \right) \vect{v}\cdot\vect{n} dw = \\
		& \int_{0}^{2\pi} R(\vect{p}(\theta)) \bar{f}\left( PDOP(\vect{p}(\theta),\vect{S}) \right) \left[ \begin{matrix}
		r_i\cos(\theta)\\
		r_i\sin(\theta)\\
		\frac{\partial r_i}{\partial s_{iz}}
		\end{matrix} \right] d\theta
	\end{aligned}
\end{equation}
上式将积分域转化为极坐标系下的相对角度，取值范围为$[0,2\pi]$。

\BiSubsection{控制策略总结与数值计算方法}{Distributed control law and numerical computation}
通过前三节的推导与分析，为了实现本章提出的主动定位服务的目的，本章所设计的节点控制策略可以直接总结为：节点需要按照式~\eqref{eq. gradienDescentUpdateRule}中的梯度下降法计算其下一时刻的最优位置，其中所用到的梯度需要按照式\eqref{eq. overviewPartialDerivative}进行展开，然后按照式~\eqref{eq. computationGradientArea} 和式~\eqref{eq. computationGradientBoundary}分别计算对应区域内的梯度值。

然而仿真应用过程发现，由于使用的节点数目相对于任务空间较为稀少，因此在机动过程中，节点在任务空间中的分布是稀疏的。因此当节点远离热点区域且其临接节点数目较少时，依靠上述过程计算得出的梯度值十分微小。此时，如果按照式~\eqref{eq. gradienDescentUpdateRule}的方法进行控制更新，节点的运动速度将近似于零，因而无法实现节点在任务空间内的有效机动。为了解决此问题，需要对式~\eqref{eq. gradienDescentUpdateRule}的更新策略进行调整。下面提出一种改进的更新策略：提取梯度的方向向量，使用固定的更新步长$\lambda$，使节点沿着梯度所指的方向运动，即：
\begin{equation} \label{eq. gradienDescentUpdateRuleDirection}
	\vect{S}_i^{k+1} = \vect{S}_i^k + \lambda \frac{\Delta \vect{S}_i^k}{|| \Delta \vect{S}_i^k ||_2}
\end{equation}
式中$\Delta \vect{S}_i^k:=\frac{\partial \mathcal{H}(\vect{S})}{\partial \vect{S}_i^k}$。上式中所做的改进，将前文计算的梯度值作为节点运动的方向，将固定的步长$\lambda$视作期望的运动速度，由此可以避免因梯度的绝对值较小而导致的问题。此外，运动速度$\lambda$的选取可以根据载体的硬件约束进行设置，从而满足实际需求。

采用式~\eqref{eq. gradienDescentUpdateRuleDirection}作为控制策略的移动协同导航服务系统，对主动定位任务的响应时间可以分为两部分。当系统获知热点区域的分布后，每个节点应当首先机动到热点区域的上空，因此称此过程所耗费的时间为机动时间。然后每当有新的节点也机动到热点区域附近时，已到达的节点群需要作出适当的队形调整以适应新的服务节点的加入。因此后一过程的需求时间也被成为调整时间。对调整时间的分析较为复杂，需要考虑各种不同初始条件和节点数目。然而机动时间的分析可以简单地通过将固定的更新步长$\lambda$视为节点的移动速度，用当前节点位置与热点区域中心点的距离除以速度来做近似地估计。

式~\eqref{eq. computationGradientArea} 和式~\eqref{eq. computationGradientBoundary}的计算涉及积分操作。在实际应用中，此类积分过程可以通过数值计算的方法获取。依照如图~\ref{fig. PDOPdiscretizationAreaBoundary}所示的离散化过程，下面将该数值计算方法简要介绍如下：
\begin{itemize}
	\item[(1)] 式~\eqref{eq. computationGradientArea}中积分域为二维平面的情况。如图~\ref{subfig. discretizationArea}所示，将任务空间离散化为$W^2$个网格，每个网格的面积为$\Delta^2$。然后使用网格中心点坐标$\vect{p}_j,1 \le j \le W^2$代表该网格内区域的坐标值。则对第$i$个节点，其覆盖区域可具体地表示为$\vect{\Omega}_i:=\{ \vect{p}_j | \forall ||\vect{p}_j - [s_{ix},s_{iy}]^T||_2 \le r_i(s_{iz}),1 \le j \le W^2 \}$。因此式~\eqref{eq. computationGradientArea}中的积分操作可以转化为：
	\begin{equation} \label{eq. computationGradientAreaNumerical}
		\begin{aligned}
			&\int_{\vect{\Omega}_i} R(\vect{p}) \frac{\partial \bar{f} \left( PDOP(\vect{p},\vect{S})\right)}{\partial \vect{S}_i} d\vect{p} = \\
			& \Delta^2 \sum_{\forall \vect{p}_j \in \vect{\Omega}_i} \left[ \begin{matrix}
			-R(\vect{p}_j) \frac{\partial \bar{f}}{\partial PDOP} tr \left[ \left( \vect{H}_j^T \vect{H}_j \right)^{-1} \frac{\partial \vect{h}_{ij}^T\vect{h}_{ij}}{\partial \vect{S}_{ix}} \left( \vect{H}_j^T \vect{H}_j \right)^{-1} \right] \\ 
			-R(\vect{p}_j) \frac{\partial \bar{f}}{\partial PDOP} tr \left[ \left( \vect{H}_j^T \vect{H}_j \right)^{-1} \frac{\partial \vect{h}_{ij}^T\vect{h}_{ij}}{\partial \vect{S}_{iy}} \left( \vect{H}_j^T \vect{H}_j \right)^{-1} \right]  \\
			-R(\vect{p}_j) \frac{\partial \bar{f}}{\partial PDOP} tr \left[ \left( \vect{H}_j^T \vect{H}_j \right)^{-1} \frac{\partial \vect{h}_{ij}^T\vect{h}_{ij}}{\partial \vect{S}_{iz}} \left( \vect{H}_j^T \vect{H}_j \right)^{-1} \right]
			\end{matrix} 
			\right]
		\end{aligned}
	\end{equation}
	式中的$\vect{H}_j$和$\vect{h}_{ij}$为式~\eqref{eq. unitDirectionVector}-\eqref{eq. collectedHMatrix}的拓展，特指在地面坐标$\vect{p}_j$处的单位方向向量和联合观测矩阵。
	
	\item[(2)] 式~\eqref{eq. computationGradientBoundary}中积分域为覆盖区域边界的情况。如图\ref{subfig. discretizationBoundary}所示，将覆盖区域的边界圆形在极坐标下离散为$U$个线段，每个线段对应夹角为$d\theta=2\pi/U$的弧,其长度为$2\pi r_i/U$。使用每个线段的中值坐标$\vect{p}_j$代表整个线段在地面上的坐标值，则边界集的定义可以被具体表示为：$\partial \vect{\Omega}_i := \{ \vect{p}_j | 1 \le j \le U \}$。此时式~\eqref{eq. computationGradientBoundary}中指标函数在覆盖边界上的梯度可以近似为：
	\begin{equation} \label{eq. computationGradientBoundaryNumerical}
		\begin{aligned}
			&\int_{\partial \vect{\Omega}_i} R(\vect{p})\bar{f}\left( PDOP(\vect{p},\vect{S}) \right) \vect{v}\cdot\vect{n} dw = \\
			& 2\pi r_i/U \sum_{\forall \vect{p}_j \in \partial \vect{\Omega}_i} R(\vect{p}_j)) \bar{f}\left( PDOP(\vect{p}_j,\vect{N}_{\vect{p}_j}^k) \right) \left[ \begin{matrix}
			r_i\cos(\theta_j)\\
			r_i\sin(\theta_j)\\
			\frac{\partial r_i}{\partial s_{iz}}
			\end{matrix} \right]
		\end{aligned}
	\end{equation}
\end{itemize}
式中$\theta_j$为坐标点$\vect{p}_j$在极坐标下的相对角度。

\begin{figure}[htbp]
	\centering	
	\subfigure{\label{subfig. discretizationArea}}\addtocounter{subfigure}{-2}
	\subfigure[Discretization over the mission space]{\subfigure[任务空间离散化]{\includegraphics[width = 0.45\textwidth]{figures_PDOP/PDOPdiscretizationCoverageArea}}}	
	\null\hfill
	\subfigure{\label{subfig. discretizationBoundary}}\addtocounter{subfigure}{-2}
	\subfigure[Discretization over the boundary]{\subfigure[覆盖区域边界离散化]{\includegraphics[width = 0.45\textwidth]{figures_PDOP/PDOPcoverageBoundary}}}
	
	\bicaption[fig. PDOPdiscretizationAreaBoundary]{}{积分的数值计算方法离散过程示意图}{Fig.$\!$}{The indication of discretization when computing the integrals}
	\vspace{-0em}
\end{figure}


\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwResult{ $\vect{S}_i^{k+1}$ }
	\KwData{ $\vect{S}_i^{k}, \vect{N}_{\vect{S}_i}^{k}, \lambda, W,U$  \;}
	
	将任务空间分割为$W^2$个网格，判定被第$i$个节点有效覆盖的网格， 将其中心点坐标收集到覆盖区域集合$\vect{\Omega}_i$中; 
	
	将覆盖边界重新分割为$U$个弧段，将所有弧段的中心点坐标收集到边界集合$\partial \vect{\Omega}_i$中;
	
	分别根据式~\eqref{eq. computationGradientAreaNumerical}和式~\eqref{eq. computationGradientBoundaryNumerical}和$\vect{N}_{\vect{S}_i}^{k}$计算指标函数在节点的覆盖区域和覆盖边界处的梯度值。
	
	根据式~\eqref{eq. overviewPartialDerivative}计算得到$\Delta \vect{S}_i^k:=\frac{\partial \mathcal{H}(\vect{S})}{\partial \vect{S}_i^k}$。
	
	
	由式~\eqref{eq. gradienDescentUpdateRuleDirection}和$\lambda$更新得到节点下一时刻的位置$\vect{S}_i^{k+1}$;
	
	\textbf{Return} $\vect{S}_i^{k+1}$
	\AlgoBiCaption{第$i$个移动节点在第$k$时刻的控制器}{Controller for $i$-th robot at timestep $k$}
	\label{alg. discontroller}
\end{algorithm}

基于以上分析，主动导航服务的控制策略可归纳为算法~\ref{alg. discontroller}。此算法的计算复杂度依赖于参数$W,U$的选取。考虑最坏的计算情况，即所有的服务节点的覆盖区域全部两两之间相互重叠，且整个任务空间被所有节点的服务范围所覆盖。此时该算法的计算复杂度正比于$O\left( 3N^2(W^2 + U) \right)$，其中$N$为系统中使用的服务节点数目。此种节点配置下，对任意地面坐标$\vect{p}$而言，其联合观测矩阵$\vect{H}$中包含所有的节点，因此复杂度中的$O\left(3N^2\right)$项对应于$\vect{H}^T\vect{H}$中的矩阵乘法。由于$\vect{H}^T\vect{H}$的矩阵大小仅为$3\times3$，因此此处的计算复杂度忽略了对该矩阵取逆的操作。

至此，本章提出的主动协同导航的分布式控制方法已经设计完毕。由算法~\ref{alg. discontroller}可看出，该算法所需要的数据全都可以在本地获得，不需要任何全局信息，因此该控制策略是分布式的，不需要中心处理节点。此外，虽然该分布式控制框架的指标函数是基于PODP分布场设计的，该过程可以十分方便地拓展到为其他精度因子，如GDOP，VDOP等。通过简单地修改，该方法也适用于传统的单纯的二维空间的覆盖问题，即仅需将式~\eqref{eq. computationGradientAreaNumerical}和式~\eqref{eq. computationGradientBoundaryNumerical}中对应于$z$轴的第三维度消去即可。

\BiSection{数值仿真与分析}{Simulation results and discussion}
本章验证前面章节提出的分布式控制算法~\ref{alg. discontroller}。后文将进行三组仿真验证实验。首先测试固定数量的节点，在不同的初始构型下，服务单个静态的热点区域。第二组实验将验证该控制策略对新节点加入和节点损失突发情况的应对处理能力。第三组仿真将验证控制策略对变化的热点区域的动态跟踪能力。

后文的三组仿真实验将统一使用如下配置。首先，导航服务的任务空间设定为$\vect{F}:=\{ [x,y,z]^T | 0\le x,y \le 100, 0 \le z \le 20 \}$，用户存在的地面空间$\vect{\Omega}$定义为$\vect{F}$中的$x-y$平面。节点对地面用户的覆盖模型选用公式~\eqref{eq. reciprocalCoverageModel}中定义的与节点飞行高度成反比例的模型，其中的参数$L$其物理意义为节点发射的导航服务信号所能传输的最远距离，以下仿真中统一设置为$L=20m$。通过修改Sigmoid函数，公式~\eqref{eq. PDOPultility1}中的光滑单调递减函数$f(\cdot)$被选取为:
\begin{equation}
	f(x) = \frac{1}{1+e^{\beta (x-\gamma)}} \quad \frac{\partial f}{\partial x} = \beta f(x)\left( f(x)-1 \right)
\end{equation}
上式中的参数$\beta$和$\gamma$被分别选定为$0.1$和$30$.式~\eqref{eq. modificationGeneralPDOP}分段函数$\bar{f}(\cdot)$中的参数设定为$a=0.05,b=0.1,c=0.15$。本节使用公式~\eqref{eq. gaussianDensity}所示的高斯密度函数表征热点区域的分布，其中参数$\alpha=1$。算法~\ref{alg. discontroller}中所需的参数为$\lambda=0.1,W=200,U=100$。

\BiSubsection{固定数量的节点服务单个静态热点}{Four pseudolites serving a static hotspot}

\begin{figure}[htbp]
	\centering
	\includegraphics[width = 0.95\textwidth]{figures_PDOP/4distribtuion1}
	\bicaption[fig. PDOP4distribution]{}{固定数量的节点服务单个静态热点区域100次的数据统计结果}{Fig.$\!$}{Statistical results for the case where four robots serve a static hotspot over 100 trials}\vspace{-0em}
\end{figure}

\begin{figure}[h]
	\centering	
	\subfigure{\label{subfig. higherPattern}}\addtocounter{subfigure}{-2}
	\subfigure[higher pattern]{\subfigure[高簇中代表构型的$\bar{f}(\cdot)$分布]{\includegraphics[width = 0.45\textwidth]{figures_PDOP/pattern1_higher}}}	
	\null\hfill
	\subfigure{\label{subfig. lowerpattern}}\addtocounter{subfigure}{-2}
	\subfigure[lower pattern]{\subfigure[低簇中代表构型的$\bar{f}(\cdot)$分布]{\includegraphics[width = 0.45\textwidth]{figures_PDOP/pattern2_lower}}}	\\
	\null\vfill
	\subfigure{\label{subfig. golbalpattern}}\addtocounter{subfigure}{-2}
	\subfigure[global pattern]{\subfigure[全局最优构型的$\bar{f}(\cdot)$分布]{\includegraphics[width = 0.45\textwidth]{figures_PDOP/globalOptimal}}}
	\bicaption[fig. PDOP4nodePatterns]{}{不同构型下PDOP场在地面热点区域附近的分布}{Fig.$\!$}{Different distribution patterns of for the $f(\cdot)$ values after reaching stable}
	\vspace{-0em}
\end{figure}

本节测试由四个移动节点组成移动协同定位服务系统，为一个位置已知的静态热点区域提供导航服务。仿真实验以空基节点从初始位置出发开始，到空基节点最终在热点区域上空形成稳态构型结束。本实验通过随机给定空基节点初始构型，通过多次打靶验证控制策略对空基集群初始构型的鲁棒性。

本节仿真实验的热点区域的密度中心点设定为$\vect{p}_d = [50,50]^T$，实验过程打靶100次，每次仿真持续1000个离散的时间步长。在每一次的仿真实验中，四个节点在任务空间中的初始位置都是随机生成的，然后使用算法~\ref{alg. discontroller}计算每个节点在下一时刻的位置。

100次打靶实验的统计结果如图~\ref{fig. PDOP4distribution}所示。该图展示了100次的打靶实验中，指标函数$\mathcal{H}(\vect{S})$随仿真时间演变的分布情况。作为对比，本节还通过采用集中式的精细化搜索算法(coarse-to-fine algorithm)，求解了四个节点在静态分布下的全局最优构型。图中顶端的水平线代表了全局最优解所能实现的最优的指标函数值。图~\ref{fig. PDOP4distribution}的子图展示了四个节点的移动导航服务系统在分布式的控制策略算法~\ref{alg. discontroller}下所能达到的最终稳定构型的分布情况。

由图中可以看出，指标函数$\mathcal{H}(\vect{S})$在全部的打靶实验中都收敛至接近全局最优的局部优化解，由此验证了本章所提控制策略的有效性，以及其对节点初始位置分布的鲁棒性。根据仿真配置，节点距离热点区域最远运动距离可以近似为$\sqrt{50^2+20^2}\approx 53.85m$，而图中显示该系统在100次的仿真实验中，其构型变化的收敛时间大约在$500$个时间步长前后，因此验证了算法~\ref{alg. discontroller}中以固定的更新步长$\lambda$作为节点运动速度的分析。

另外通过观察子图可以看出，100次的仿真实验最终稳定地分布在不同的局部最优构型上。该现象的产生是由于本章设计的控制器基于梯度下的更新策略，因此其仅能达到局部最优的系统构型，而无法保证得到全局最优解。100次打靶实验的稳态构型主要集中于两个分布簇中，即处于$5.247-5.248$区域的高簇和处于$5.245-5.247$区域的低簇。这两个分布簇中的代表构型以及全局最优构型在地面任务空间的PDOP场如图~\ref{fig. PDOP4nodePatterns}所示，图中的颜色深浅代表在该地面坐标上$\bar{f}(PODP(\vect{S},\vect{p}))$的数值,该数值越大，颜色越明亮，该位置获得的PDOP值越小，代表定位信号的构型质量越好。通过对比其PDOP场的分布可知，局部最优解往往导致热点区域内的一部分地面坐标与四个节点的分布接近共面，从而使得其整体的表现不如全局最优解。

%最后，图~\ref{fig. PDOP4distribution}显示，全部100次的仿真的性能表现全部都是局部最优解，即没有一次仿真最终稳定地收敛到全局最优解。然而，通过第三个实验仿真可以验证，如果节点的初始位置选取合适，本章所提出的分布式控制策略是可以使服务系统达到到全局最优解。

\BiSubsection{新节点加入与节点损失}{The loss and resurgence of pseudolites}
本节验证设计的分布式控制策略应对新节点加入和节点损失等特殊情况的控制性能。该实验的仿真流程设计如下：
\begin{itemize}
	\item[(1)] 第一阶段：仿真开始时，随机初始化三个节点，使用设计的分布式控制策略，为单个已知的静态热点区域提供定位导航服务。此过程持续700个时间步长。
	
	\item[(2)] 第二阶段： 当三个节点达到稳定的构型后，将两个新节点随机地放置于任务空间之中来表征新服务节点的加入。等待700个时间步长，让所有5个节点重新达到稳态。
	
	\item[(3)] 第三阶段：最后从5个节点中随机去掉一个节点，以模拟在特殊情况下节点损失情况，此过程依然持续700步。
\end{itemize}

上述过程中存在三处随机事件，即第一阶段三个节点的初始位置，第二阶段新加入的两个节点的初始位置，以及第三阶段损失节点的身份。为了验证上述随机事件对构型控制方法~\ref{alg. discontroller}的影响，本实验进行30次打靶。在每次的打靶实验中，上述三处事件都随机产生。30次打靶实验得到的指标函数$\mathcal{H}(\vect{S})$随时间变化的分布如图~\ref{fig. Nodefailurerecurren}所示。图~\ref{fig. Nodefailurerecurren}中的子图分别显示指标函数在不同时刻的变化细节。其中子图(A)显示了指标函数在排除一组性能最坏的仿真后的具体分布;子图(B)展示了当两个新节点加入后，指标函数的增加变化过程;子图(C)描述了当排除了节点损失导致的性能最坏的仿真组后，指标函数的具体波动;子图(D)给出了节点数目变化后的最终稳态构型分布。

\begin{figure}[htbp]
	\centering
	\includegraphics[width = 0.95\textwidth]{figures_PDOP/Nodefailurerecurren1}
	\bicaption[fig. Nodefailurerecurren]{}{30次节点增加和减少仿真对拒止区域的PDOP场覆盖效果的数据统计结果}{Fig.$\!$}{Statistical results for the case node loss and resurgence over 30 trials}\vspace{-0.5em}
\end{figure}

第一阶段的仿真实际上是前一节仿真结果在三个服务节点系统中的拓展。与前一节的结果有所不同的是，三个节点组成的服务系统出现了一组奇异解，即图~\ref{fig. Nodefailurerecurren} (A)中虚线所示的情况。此时，三个节点在稳定状态下收敛到一个三点共线的奇异构型，从而导致地面的热点区域无法接收到有效的导航服务。排除这一特殊仿真结果后，多数的仿真结果可以稳定到一个有效的局部最优构型。

图~\ref{fig. Nodefailurerecurren} (B)显示，新节点的加入从整体上提升了系统的导航服务能力。从图中可以观察得到，5个节点的服务系统在达到稳定构型后，其对应的指标函数的分布比3个节点系统的稳态分布更为紧凑，整体数值更高。即使是子图(A)中的奇异实验组，在新节点加入后也迅速恢复了有效的导航服务。

子图(C)显示，在一个节点被随机去除后，系统在节点损失的瞬间产生了覆盖指标下降的波动，因此节点损失会降低导航服务系统的服务质量。子图(C)中的虚线展示了全部30次仿真之中对节点损失表现最差的一组实验，即节点损失后剩余节点恰好处在共线的奇异配置下。然而从图中可以看出，剩下的四个节点在节点损失之后立刻在分布式控制策略的机动下脱离奇异配置，从新形成了一个有效的导航服务构型。

最后从子图(D)中可以看出，经过节点损失而形成的四个节点的服务系统，其稳态构型最终也趋向于两个分布簇。子图(D)中高分布簇的位置与前一节仿真结果中高分布簇的位置近似相同，然而低分布簇的位置却比前一节低分布簇的位置更低。由此可以总结出，经过节点损失后的系统，会比相同配置的系统从初始位置进行机动的情况获得更差的性能表现。

该仿真实验同时也验证了本章所提出的分布式控制策略对新节点的加入以及节点损失的鲁棒性。残存节点对新节点加入和节点损失这样的特殊事件的反映并不强烈。这一现象的背后是梯度下降法的局部最优特性。由于多节点系统$N\ge4$存在许多局部最优解，因此当特殊事件发生后，残存的节点系统只需要迁徙到其距离当前状态最近的局部解即可，因此无需大规模内部构型的调整。

\BiSubsection{对动态切换热点的跟踪服务}{The case for switching hotspots}

最后一组实验验证控制策略对动态变化的热点区域的追踪服务能力。该实验的仿真流程如下：

\begin{itemize}
	\item[(1)] 初始条件下三个节点随机分布于任务空间中。
	\item[(2)] 在整个仿真过程中，热点区域的密度中心点不再是固定的。此次仿真实验一共设定有四个不同的密度中心，分别为$d_{c1} = [20,20]^T,d_{c2} = [80,20]^T,d_{c3} = [80,80]^T,d_{c4} = [20,80]^T$。每隔500个仿真间隔，热点区域的中心点将切换至下一个。
\end{itemize}

\begin{figure}[htbp]
	\centering	
	\subfigure{\label{subfig. PDOPswitchedDistribution}}\addtocounter{subfigure}{-2}
	\subfigure[Joint distribtuion of PDOP filed after reaching stable for all four hotspots]{\subfigure[四个热点处稳态时的PDOP场分布情况]{\includegraphics[width = 0.48\textwidth]{figures_PDOP/PDOPdistributionFourPoints}}}	
	\null\hfill
	\subfigure{\label{subfig. PDOPswitchedTrajectories}}\addtocounter{subfigure}{-2}
	\subfigure[Trajectories for all four agents when switching among hotspots]{\subfigure[四个智能体的转移路径]{\includegraphics[width = 0.48\textwidth]{figures_PDOP/Trajectories3DFourPoints}}}
	
	\bicaption[fig. PDOP3node4hotspots]{}{切换热点的构型控制方法效果}{Fig.$\!$}{The performance of the configuration control serving switching hotspots}
	\vspace{-0.5em}
\end{figure}

与前两个实验有所不同的是，因为热点的切换间隔仅为500个时间步长，因此需要设置算法~\ref{alg. discontroller}中固定更新步长为$\lambda=0.2$以加快节点的运动速度，使得节点在热点切换后，能够及时地从前一个热点区域上空机动到新的热点区域。一组典型的实验仿真结果如图~\ref{fig. PDOP3node4hotspots}所示。

图~\ref{subfig. PDOPswitchedDistribution}是融合了三个节点的导航系统在不同的稳态构型下对四个个热点区域的PDOP场服务分布。其中左下角对第一个热点区域的PDOP场分布与全局最优解所能提供的PDOP场分布完全相同。因此说明了本章所提出的分布式控制策略在合适的初始配置下能够收敛于全局最优解。之后的三个分布全部退化为局部最优解。图~\ref{subfig. PDOPswitchedTrajectories}展示了所有节点对四个热点区域服务的切换轨迹。除了在第一个热点区域上空的轨迹变化较为曲折外，向其他三个热点的转移过程以及节点抵达后节点内部构型的调整过程全都较为平滑。

仿真结果验证了分布式控制策略~\ref{alg. discontroller}可以实现对动态变化的热点区域进行跟踪，并具备自动地形成基于PDOP场分布的最优服务构型的能力。

\BiSection{小结}{Summary}

本章研究了移动协同定位系统的主动优化导航服务的问题，提出了一种基于梯度下降法的分布式构型控制方法，使得移动协同定位系统具备自动机动到指定热点区域上空，自主改变内部构型从而向热点区域内的地面用户提供优化定位服务的能力。本章通过严格的理论推导，给出了以PDOP分布场为量化指标的导航服务质量对节点位置变化的梯度计算公式。分别给出了梯度在节点有效覆盖范围内和覆盖边界处的具体计算方法。由于本章设计的控制策略是分布式的，理论分析和仿真实验结果表明，使用该控制策略，移动协同导航系统除了能够提供主动导航服务外，还具备如下特性：

\begin{itemize}
	\item[(1)] 对节点的初始位置分布具有鲁棒性。
	\item[(2)] 无中心计算节点，因此抗节点损失，可应用于不同规模的系统。
	\item[(3)] 具备对热点区域的动态追踪能力。
	\item[(4)] 可以通过调整更新步长改变节点的运动速度以及系统的收敛时间。
\end{itemize}

本章提出的控制策略的局限性有：
\begin{itemize}
	\item[(1)] 此控制策略需要节点获取自身的精准定位信息。然而通过协同定位或者SLAM等方法获取的定位信息总是不精确，存在定位误差。
	\item[(2)] 当前的控制策略假设热点区域的位置和分布对任一节点在任何时刻都是已知的。
	\item[(3)] 本章的优化方法是基于梯度下降法进行求解的，一般只能获取局部最优解。
\end{itemize}







